using DataFrames
using QuadGK
using Distributions
using JLD
using Optim
using ForwardDiff
using CSV
using LinearAlgebra
using SpecialFunctions
using Random
using DelimitedFiles
using SparseArrays
using Interpolations

clearconsole()

#-------------------------------------------------------------
# include Functions
#-------------------------------------------------------------
#cd("$(pwd())/Dropbox/Heterogeneity/Software/KS_Simulation/")
readDir = "$(pwd())/Functions/"
include(readDir *"vech.jl");
include(readDir *"logSpline_Procedures.jl");
include(readDir *"VAR_Procedures.jl");
include(readDir *"Loaddata.jl");
include(readDir *"EmpPercentiles_Procedures.jl")

#-------------------------------------------------------------
# choose specification files
#-------------------------------------------------------------
nfVARSpec = "3"
Tend      = 260
K         = 8

specDir   = "$(pwd())/SpecFiles/"
include(specDir * "/fVARspec" * nfVARSpec * ".jl")
nKSpec    = "K$(K)_"

#-------------------------------------------------------------
# load aggregate data (Zt and Kexact)
#-------------------------------------------------------------
nDataSel = "2"
dataDir  = "$(pwd())/data/Para" * nDataSel *"/"
agg_data    = loadaggdata(1,Tend)
period_agg  = size(agg_data)[1]
n_agg       = size(agg_data)[2]

#-------------------------------------------------------------
# Load coefficients from density estimation
#-------------------------------------------------------------
sName   = "fVAR" * nfVARSpec
loaddir  = "$(pwd())/results/Para" *nDataSel* "/" * sName *"/";

knots_all    = readdlm(loaddir * sName * "_knots_all.csv", ',', Float64, '\n'; skipstart=1);
PhatDensCoef = readdlm(loaddir * nKSpec * sName * "_PhatDensCoef.csv", ',', Float64, '\n'; skipstart=1);
PhatDensCoef = PhatDensCoef[1:Tend,:]

#-------------------------------------------------------------
# Select Knots
#-------------------------------------------------------------

ii=getindex.(findall(K_vec.-K.==0),[1 2])[1] # find index ii where K==K_vec
knots = knots_all[quant_sel[ii,:].==1]

#-------------------------------------------------------------
# Initializations
#-------------------------------------------------------------

# compute empirical CDF for a grid of values
ngrid     = 20
grid_temp = range(0.2, stop=2, length=ngrid);
# ngrid     = 50
# grid_temp = range(0.05, stop=2.5, length=ngrid);
grid_temp = [0.0; grid_temp]

# initalize the matrices for the time series of percentiles
vec_percs = [0.1; 0.2; 0.5; 0.8; 0.9]
emp_percs_all = zeros(Tend,length(vec_percs))

#-------------------------------------------------------------
# Loop over time periods to compute percentiles
#-------------------------------------------------------------

for tt = 1:Tend

    print("CDF and Percentiles for Period = $tt \n")
    unrate    = 0.0
    @time emp_percs = DensPercentiles(PhatDensCoef[tt,:]', knots, unrate, xgrid, grid_temp, vec_percs)
    emp_percs_all[tt,:] = emp_percs'

end

savedir  = "$(pwd())/results/Para" *nDataSel* "/" * sName *"/";
CSV.write(savedir * nKSpec * sName * "_PredPctL_MLE.csv", DataFrame(emp_percs_all,:auto));
